Figure 2: Core Microbiome Analysis

Recommendation for reading and cutoffs: https://www.pnas.org/doi/10.1073/pnas.2104429118

Overview: 1. PCoA Analysis - complete microbiome using Unifrac? distances 2. Alpha diversity and betadispersion 3. Determine the core across control samples 4. How much of the community is captured? 5. How does community change with exposure to ABS? 6. Which core members are differential abundant at TP4,TP5,TP6?

library(tidyr)
library(tidyverse)

── Attaching core tidyverse packages ──────────────────────────────────────────────────────────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.3     ✔ purrr     1.0.2
✔ forcats   1.0.0     ✔ readr     2.1.4
✔ ggplot2   3.4.4     ✔ stringr   1.5.0
✔ lubridate 1.9.3     ✔ tibble    3.2.1── Conflicts ────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors

library(ggplot2)
library(vegan)

Loading required package: permute
Loading required package: lattice
This is vegan 2.6-4

library(phyloseq)
library(ggrepel)
library(patchwork)
library(microbiome)

microbiome R package (microbiome.github.com)
    


 Copyright (C) 2011-2022 Leo Lahti, 
    Sudarshan Shetty et al. <microbiome.github.io>


Attaching package: ‘microbiome’

The following object is masked from ‘package:vegan’:

    diversity

The following object is masked from ‘package:ggplot2’:

    alpha

The following object is masked from ‘package:base’:

    transform

library(RColorBrewer)

plotting conditions

trt16<-c('ABS1','ABS2','Control','Experimental', 'FASW','Priming','blank', 'FAASW' )
l2p16<-c('indianred1','deeppink4','gray16','burlywood','burlywood','burlywood', 'gray81', 'darkgreen')
names(l2p16) = trt16

shpsc<-c(0,1,3,4,8,15,16)
names(shpsc)=c('T0','T1','T2','T3','T4','T5','T6')

theme_new<- theme_bw() + theme(panel.grid.major = element_line(color='gray90'), panel.grid.minor = element_line(color='gray90'), text= element_text(size = 15), axis.text.x = element_text(angle = 65,hjust = 1), legend.text = element_text(size = 20), axis.title = element_text(size = 15), title = element_text(size = 20))

read in cleaned data

ps<-readRDS(file = 'ps.cleantree.noblanks.rds')
ps<-prune_samples(!(sample_names(ps) %in% c('EB1_S94', 'EB1_S95','EB2_S95')),ps) #remove EB samples 
ps

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 414 taxa and 82 samples ]
sample_data() Sample Data:       [ 82 samples by 8 sample variables ]
tax_table()   Taxonomy Table:    [ 414 taxa by 7 taxonomic ranks ]
phy_tree()    Phylogenetic Tree: [ 414 tips and 412 internal nodes ]
refseq()      DNAStringSet:      [ 414 reference sequences ]

Fig 2A: Alpha Diveristy

ps.red<-subset_samples(ps, Time %in% c('T0','T1','T4','T5','T6'))

p<-plot_richness(ps.red, measures=c("Observed"))
p_df<-p$data

a<-ggplot(p_df, aes(y=value, x=Condition, fill=Condition)) + 
  geom_boxplot(alpha=0.5, outlier.shape = NA, color='black') + 
   scale_fill_manual(values=l2p16)+scale_color_manual(values=l2p16)+
  geom_jitter(width=0.1,height=0,size=1, alpha=0.5) + facet_grid(~Time,scales='free',space='free')+ 
  theme_new + ylab('Nr. ASVs') + xlab('') + guides(fill=F, group=F,color=F) + ylim(c(0,140)) + 
  theme(axis.text.x = element_blank(),panel.grid = element_blank(),panel.grid.minor = element_blank(),rect  = element_blank())
a

stats

p_df.trt<-p_df 
hist(p_df.trt$value)
p_df.trt$Time<-as.factor(p_df.trt$Time)
p_df.trt$Condition<-as.factor(p_df.trt$Condition)


aov.trt<-aov(value ~ Condition*Time, p_df.trt)
shapiro.test(residuals(aov.trt))

    Shapiro-Wilk normality test

data:  residuals(aov.trt)
W = 0.99104, p-value = 0.9574

summary(aov.trt)

               Df Sum Sq Mean Sq F value   Pr(>F)    
Condition       3   7147  2382.2   8.494 0.000159 ***
Time            4   6159  1539.6   5.490 0.001199 ** 
Condition:Time  4   9508  2377.0   8.475 4.16e-05 ***
Residuals      42  11780   280.5                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

aov.trt

Call:
   aov(formula = value ~ Condition * Time, data = p_df.trt)

Terms:
                Condition      Time Condition:Time Residuals
Sum of Squares   7146.709  6158.544       9507.906 11779.600
Deg. of Freedom         3         4              4        42

Residual standard error: 16.74714
8 out of 20 effects not estimable
Estimated effects may be unbalanced

par(mfrow=(c(2,2)))

plot(aov.trt)

data_summary <- group_by(p_df.trt, Condition, Time) %>%
  summarise(mean=mean(value), sd=sd(value), se=sd(value)/sqrt(length(value))) %>%
  arrange(desc(mean))

`summarise()` has grouped output by 'Condition'. You can override using the `.groups` argument.

print(data_summary)

tukey <- TukeyHSD(aov.trt)
print(tukey)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = value ~ Condition * Time, data = p_df.trt)

$Condition
                          diff        lwr      upr     p adj
ABS2-ABS1             8.794872  -8.180495 25.77024 0.5150285
Control-ABS1         27.747253  11.937845 43.55666 0.0001624
Experimental-ABS1    21.861538  -1.712666 45.43574 0.0777548
Control-ABS2         18.952381   3.807933 34.09683 0.0090201
Experimental-ABS2    13.066667 -10.066859 36.20019 0.4401964
Experimental-Control -5.885714 -28.177738 16.40631 0.8940480

$Time
            diff       lwr        upr     p adj
T1-T0  -5.841270 -37.65865 25.9761110 0.9845071
T4-T0  -3.384334 -33.95348 27.1848171 0.9977642
T5-T0  -9.107274 -39.47103 21.2564829 0.9115463
T6-T0 -27.233455 -57.41807  2.9511624 0.0943868
T4-T1   2.456936 -18.23848 23.1523491 0.9970659
T5-T1  -3.266004 -23.65681 17.1247993 0.9907266
T6-T1 -21.392186 -41.51526 -1.2691070 0.0322421
T5-T4  -5.722940 -24.10533 12.6594456 0.9000293
T6-T4 -23.849122 -41.93408 -5.7641683 0.0045120
T6-T5 -18.126182 -35.86175 -0.3906164 0.0430011

$`Condition:Time`
                                 diff         lwr        upr     p adj
ABS2:T0-ABS1:T0                    NA          NA         NA        NA
Control:T0-ABS1:T0                 NA          NA         NA        NA
Experimental:T0-ABS1:T0            NA          NA         NA        NA
ABS1:T1-ABS1:T0                    NA          NA         NA        NA
ABS2:T1-ABS1:T0                    NA          NA         NA        NA
Control:T1-ABS1:T0                 NA          NA         NA        NA
Experimental:T1-ABS1:T0            NA          NA         NA        NA
ABS1:T4-ABS1:T0                    NA          NA         NA        NA
ABS2:T4-ABS1:T0                    NA          NA         NA        NA
Control:T4-ABS1:T0                 NA          NA         NA        NA
Experimental:T4-ABS1:T0            NA          NA         NA        NA
ABS1:T5-ABS1:T0                    NA          NA         NA        NA
ABS2:T5-ABS1:T0                    NA          NA         NA        NA
Control:T5-ABS1:T0                 NA          NA         NA        NA
Experimental:T5-ABS1:T0            NA          NA         NA        NA
ABS1:T6-ABS1:T0                    NA          NA         NA        NA
ABS2:T6-ABS1:T0                    NA          NA         NA        NA
Control:T6-ABS1:T0                 NA          NA         NA        NA
Experimental:T6-ABS1:T0            NA          NA         NA        NA
Control:T0-ABS2:T0                 NA          NA         NA        NA
Experimental:T0-ABS2:T0            NA          NA         NA        NA
ABS1:T1-ABS2:T0                    NA          NA         NA        NA
ABS2:T1-ABS2:T0                    NA          NA         NA        NA
Control:T1-ABS2:T0                 NA          NA         NA        NA
Experimental:T1-ABS2:T0            NA          NA         NA        NA
ABS1:T4-ABS2:T0                    NA          NA         NA        NA
ABS2:T4-ABS2:T0                    NA          NA         NA        NA
Control:T4-ABS2:T0                 NA          NA         NA        NA
Experimental:T4-ABS2:T0            NA          NA         NA        NA
ABS1:T5-ABS2:T0                    NA          NA         NA        NA
ABS2:T5-ABS2:T0                    NA          NA         NA        NA
Control:T5-ABS2:T0                 NA          NA         NA        NA
Experimental:T5-ABS2:T0            NA          NA         NA        NA
ABS1:T6-ABS2:T0                    NA          NA         NA        NA
ABS2:T6-ABS2:T0                    NA          NA         NA        NA
Control:T6-ABS2:T0                 NA          NA         NA        NA
Experimental:T6-ABS2:T0            NA          NA         NA        NA
Experimental:T0-Control:T0         NA          NA         NA        NA
ABS1:T1-Control:T0                 NA          NA         NA        NA
ABS2:T1-Control:T0                 NA          NA         NA        NA
Control:T1-Control:T0             1.5  -46.805061  49.805061 1.0000000
Experimental:T1-Control:T0      -17.6  -63.788473  28.588473 0.9946379
ABS1:T4-Control:T0                1.0  -50.640283  52.640283 1.0000000
ABS2:T4-Control:T0              -16.8  -62.988473  29.388473 0.9968728
Control:T4-Control:T0           -28.2  -74.388473  17.988473 0.7202862
tkl         NA          NA         NA        NA
ABS1:T5-Control:T0              -51.4  -97.588473  -5.211527 0.0160516
ABS2:T5-Control:T0              -14.8  -60.988473  31.388473 0.9993584
Control:T5-Control:T0            -7.5  -55.805061  40.805061 1.0000000
Experimental:T5-Control:T0         NA          NA         NA        NA
ABS1:T6-Control:T0              -51.8  -97.988473  -5.611527 0.0146500
ABS2:T6-Control:T0              -60.4 -106.588473 -14.211527 0.0018482
Control:T6-Control:T0           -16.2  -62.388473  29.988473 0.9979815
Experimental:T6-Control:T0         NA          NA         NA        NA
ABS1:T1-Experimental:T0            NA          NA         NA        NA
ABS2:T1-Experimental:T0            NA          NA         NA        NA
Control:T1-Experimental:T0         NA          NA         NA        NA
Experimental:T1-Experimental:T0    NA          NA         NA        NA
ABS1:T4-Experimental:T0            NA          NA         NA        NA
ABS2:T4-Experimental:T0            NA          NA         NA        NA
Control:T4-Experimental:T0         NA          NA         NA        NA
Experimental:T4-Experimental:T0    NA          NA         NA        NA
ABS1:T5-Experimental:T0            NA          NA         NA        NA
ABS2:T5-Experimental:T0            NA          NA         NA        NA
Control:T5-Experimental:T0         NA          NA         NA        NA
Experimental:T5-Experimental:T0    NA          NA         NA        NA
ABS1:T6-Experimental:T0            NA          NA         NA        NA
ABS2:T6-Experimental:T0            NA          NA         NA        NA
Control:T6-Experimental:T0         NA          NA         NA        NA
Experimental:T6-Experimental:T0    NA          NA         NA        NA
ABS2:T1-ABS1:T1                    NA          NA         NA        NA
Control:T1-ABS1:T1                 NA          NA         NA        NA
Experimental:T1-ABS1:T1            NA          NA         NA        NA
ABS1:T4-ABS1:T1                    NA          NA         NA        NA
ABS2:T4-ABS1:T1                    NA          NA         NA        NA
Control:T4-ABS1:T1                 NA          NA         NA        NA
Experimental:T4-ABS1:T1            NA          NA         NA        NA
ABS1:T5-ABS1:T1                    NA          NA         NA        NA
ABS2:T5-ABS1:T1                    NA          NA         NA        NA
Control:T5-ABS1:T1                 NA          NA         NA        NA
Experimental:T5-ABS1:T1            NA          NA         NA        NA
ABS1:T6-ABS1:T1                    NA          NA         NA        NA
ABS2:T6-ABS1:T1                    NA          NA         NA        NA
Control:T6-ABS1:T1                 NA          NA         NA        NA
Experimental:T6-ABS1:T1            NA          NA         NA        NA
Control:T1-ABS2:T1                 NA          NA         NA        NA
Experimental:T1-ABS2:T1            NA          NA         NA        NA
ABS1:T4-ABS2:T1                    NA          NA         NA        NA
ABS2:T4-ABS2:T1                    NA          NA         NA        NA
Control:T4-ABS2:T1                 NA          NA         NA        NA
Experimental:T4-ABS2:T1            NA          NA         NA        NA
ABS1:T5-ABS2:T1                    NA          NA         NA        NA
ABS2:T5-ABS2:T1                    NA          NA         NA        NA
Control:T5-ABS2:T1                 NA          NA         NA        NA
Experimental:T5-ABS2:T1            NA          NA         NA        NA
ABS1:T6-ABS2:T1                    NA          NA         NA        NA
ABS2:T6-ABS2:T1                    NA          NA         NA        NA
Control:T6-ABS2:T1                 NA          NA         NA        NA
Experimental:T6-ABS2:T1            NA          NA         NA        NA
Experimental:T1-Control:T1      -19.1  -61.526821  23.326821 0.9698630
ABS1:T4-Control:T1               -0.5  -48.805061  47.805061 1.0000000
ABS2:T4-Control:T1              -18.3  -60.726821  24.126821 0.9800159
Control:T4-Control:T1           -29.7  -72.126821  12.726821 0.4918565
Experimental:T4-Control:T1         NA          NA         NA        NA
ABS1:T5-Control:T1              -52.9  -95.326821 -10.473179 0.0037057
ABS2:T5-Control:T1              -16.3  -58.726821  26.126821 0.9941168
Control:T5-Control:T1            -9.0  -53.721797  35.721797 0.9999994
Experimental:T5-Control:T1         NA          NA         NA        NA
ABS1:T6-Control:T1              -53.3  -95.726821 -10.873179 0.0033305
ABS2:T6-Control:T1              -61.9 -104.326821 -19.473179 0.0003099
Control:T6-Control:T1           -17.7  -60.126821  24.726821 0.9857164
Experimental:T6-Control:T1         NA          NA         NA        NA
ABS1:T4-Experimental:T1          18.6  -27.588473  64.788473 0.9901418
ABS2:T4-Experimental:T1           0.8  -39.200391  40.800391 1.0000000
Control:T4-Experimental:T1      -10.6  -50.600391  29.400391 0.9999546
Experimental:T4-Experimental:T1    NA          NA         NA        NA
ABS1:T5-Experimental:T1         -33.8  -73.800391   6.200391 0.1922078
ABS2:T5-Experimental:T1           2.8  -37.200391  42.800391 1.0000000
Control:T5-Experimental:T1       10.1  -32.326821  52.526821 0.9999910
Experimental:T5-Experimental:T1    NA          NA         NA        NA
ABS1:T6-Experimental:T1         -34.2  -74.200391   5.800391 0.1778735
ABS2:T6-Experimental:T1         -42.8  -82.800391  -2.799609 0.0250200
Control:T6-Experimental:T1        1.4  -38.600391  41.400391 1.0000000
Experimental:T6-Experimental:T1    NA          NA         NA        NA
ABS2:T4-ABS1:T4                 -17.8  -63.988473  28.388473 0.9939096
Control:T4-ABS1:T4              -29.2  -75.388473  16.988473 0.6667489
Experimental:T4-ABS1:T4            NA          NA         NA        NA
ABS1:T5-ABS1:T4                 -52.4  -98.588473  -6.211527 0.0127614
ABS2:T5-ABS1:T4                 -15.8  -61.988473  30.388473 0.9985179
Control:T5-ABS1:T4               -8.5  -56.805061  39.805061 0.9999999
Experimental:T5-ABS1:T4            NA          NA         NA        NA
ABS1:T6-ABS1:T4                 -52.8  -98.988473  -6.611527 0.0116322
ABS2:T6-ABS1:T4                 -61.4 -107.588473 -15.211527 0.0014377
Control:T6-ABS1:T4              -17.2  -63.388473  28.988473 0.9958800
Experimental:T6-ABS1:T4            NA          NA         NA        NA
Control:T4-ABS2:T4              -11.4  -51.400391  28.600391 0.9998702
Experimental:T4-ABS2:T4            NA          NA         NA        NA
ABS1:T5-ABS2:T4                 -34.6  -74.600391   5.400391 0.1643710
ABS2:T5-ABS2:T4                   2.0  -38.000391  42.000391 1.0000000
Control:T5-ABS2:T4                9.3  -33.126821  51.726821 0.9999975
Experimental:T5-ABS2:T4            NA          NA         NA        NA
ABS1:T6-ABS2:T4                 -35.0  -75.000391   5.000391 0.1516798
ABS2:T6-ABS2:T4                 -43.6  -83.600391  -3.599609 0.0203729
Control:T6-ABS2:T4                0.6  -39.400391  40.600391 1.0000000
Experimental:T6-ABS2:T4            NA          NA         NA        NA
Experimental:T4-Control:T4         NA          NA         NA        NA
ABS1:T5-Control:T4              -23.2  -63.200391  16.800391 0.7901962
ABS2:T5-Control:T4               13.4  -26.600391  53.400391 0.9988609
Control:T5-Control:T4            20.7  -21.726821  63.126821 0.9387702
Experimental:T5-Control:T4         NA          NA         NA        NA
ABS1:T6-Control:T4              -23.6  -63.600391  16.400391 0.7681937
ABS2:T6-Control:T4              -32.2  -72.200391   7.800391 0.2581287
Control:T6-Control:T4            12.0  -28.000391  52.000391 0.9997349
Experimental:T6-Control:T4         NA          NA         NA        NA
ABS1:T5-Experimental:T4            NA          NA         NA        NA
ABS2:T5-Experimental:T4            NA          NA         NA        NA
Control:T5-Experimental:T4         NA          NA         NA        NA
Experimental:T5-Experimental:T4    NA          NA         NA        NA
ABS1:T6-Experimental:T4            NA          NA         NA        NA
ABS2:T6-Experimental:T4            NA          NA         NA        NA
Control:T6-Experimental:T4         NA          NA         NA        NA
Experimental:T6-Experimental:T4    NA          NA         NA        NA
ABS2:T5-ABS1:T5                  36.6   -3.400391  76.600391 0.1085233
Control:T5-ABS1:T5               43.9    1.473179  86.326821 0.0356593
Experimental:T5-ABS1:T5            NA          NA         NA        NA
ABS1:T6-ABS1:T5                  -0.4  -40.400391  39.600391 1.0000000
ABS2:T6-ABS1:T5                  -9.0  -49.000391  31.000391 0.9999963
Control:T6-ABS1:T5               35.2   -4.800391  75.200391 0.1456312
Experimental:T6-ABS1:T5            NA          NA         NA        NA
Control:T5-ABS2:T5                7.3  -35.126821  49.726821 1.0000000
Experimental:T5-ABS2:T5            NA          NA         NA        NA
ABS1:T6-ABS2:T5                 -37.0  -77.000391   3.000391 0.0994925
ABS2:T6-ABS2:T5                 -45.6  -85.600391  -5.599609 0.0120320
Control:T6-ABS2:T5               -1.4  -41.400391  38.600391 1.0000000
Experimental:T6-ABS2:T5            NA          NA         NA        NA
Experimental:T5-Control:T5         NA          NA         NA        NA
ABS1:T6-Control:T5              -44.3  -86.726821  -1.873179 0.0324703
ABS2:T6-Control:T5              -52.9  -95.326821 -10.473179 0.0037057
Control:T6-Control:T5            -8.7  -51.126821  33.726821 0.9999992
Experimental:T6-Control:T5         NA          NA         NA        NA
ABS1:T6-Experimental:T5            NA          NA         NA        NA
ABS2:T6-Experimental:T5            NA          NA         NA        NA
Control:T6-Experimental:T5         NA          NA         NA        NA
Experimental:T6-Experimental:T5    NA          NA         NA        NA
ABS2:T6-ABS1:T6                  -8.6  -48.600391  31.400391 0.9999982
Control:T6-ABS1:T6               35.6   -4.400391  75.600391 0.1341122
Experimental:T6-ABS1:T6            NA          NA         NA        NA
Control:T6-ABS2:T6               44.2    4.199609  84.200391 0.0174282
Experimental:T6-ABS2:T6            NA          NA         NA        NA
Experimental:T6-Control:T6         NA          NA         NA        NA

head(p_df)

p_df %>% group_by(Time, Condition) %>% summarize(mean_asvs=mean(value), se_asv=sd(value)/sqrt(length(value)))

`summarise()` has grouped output by 'Time'. You can override using the `.groups` argument.

#replace NAs with 0s to get multcomp to run 
tukey$`Condition:Time`[!complete.cases(tukey$`Condition:Time`),] <- 0

library(multcompView)
tukey.cld <-multcompLetters4(aov.trt, tukey)
print(tukey.cld)

$Condition
     Control Experimental         ABS2         ABS1 
         "a"         "ab"          "b"          "b" 

$Time
  T0   T1   T4   T5   T6 
"ab"  "a"  "a"  "a"  "b" 

$`Condition:Time`
     Control:T1         ABS1:T4      Control:T0      Control:T5         ABS2:T5      Control:T6         ABS2:T4 Experimental:T1 
            "a"             "a"             "a"             "a"            "ab"            "ab"            "ab"            "ab" 
     Control:T4         ABS1:T5         ABS1:T6         ABS2:T6         ABS1:T0         ABS2:T0 tkl         ABS1:T1 
          "abc"            "bc"            "bc"             "c"             "d"             "e"             "f"             "g" 
        ABS2:T1 Experimental:T4 Experimental:T5 Experimental:T6 
            "h"             "i"             "j"             "k" 

tkl<-as.vector(tukey.cld$`Condition:Time`)
tkl.df<-data.frame(trt=tkl['Letters'])
tkl.df<-tkl.df %>% rownames_to_column() %>% separate(rowname,into=c('Condition','Time'))
head(tkl.df)

tkl.df$Condition<-factor(tkl.df$Condition, ordered=T, levels=c('Control','Experimental','ABS1','ABS2'))
tkl.df<-tkl.df[!tkl.df$Letters %in% c('d','e','f','g','h','i','j','k'),]

head(p_df)

p_df$Condition<-factor(p_df$Condition, ordered=T, levels=c('Control','Experimental','ABS1','ABS2'))

a<-ggplot(p_df, aes(y=value, x=Condition, fill=Condition)) + 
  geom_boxplot(alpha=0.5, outlier.shape = NA, color='black') + 
  geom_text(data=tkl.df, aes(y=10,x=Condition,label=Letters))+
   scale_fill_manual(values=l2p16)+scale_color_manual(values=l2p16)+
  geom_jitter(width=0.1,height=0,size=1, alpha=0.5) + facet_grid(~Time,scales='free',space='free')+ 
  theme_new + ylab('Nr. ASVs') + xlab('') + guides(fill=F, group=F,color=F) + ylim(c(0,140)) + 
  theme(axis.text.x = element_blank(),panel.grid = element_blank(),panel.grid.minor = element_blank(),rect  = element_blank())
a

Indentify the Core Analyisis

ASVs present in at least 90% of control samples at any relative abundance

# subset for control group
ps.control<-subset_samples(ps.red, Condition %in% c('Control'))

#removes taxa with 0 counts
ps.control <- prune_taxa(taxa_sums(ps.control) > 0, ps.control)

#transforms to rel. abundance 
pseq.rel <- microbiome::transform(ps.control, "compositional")

plots across prevalance thresholds to help guide selection of the core

# Core with compositionals:
prevalences <- seq(.05, 1, .05)
detections <- round(10^seq(log10(1e-3), log10(.2), length = 10), 3)

# Also define gray color palette
gray <- gray(seq(0,1,length=5))
# Only include core taxa
pseq.core <- core(pseq.rel, detection =0, prevalence = 90/100)
p <- plot_core(pseq.core, plot.type = "heatmap",
           colours = rev(brewer.pal(5, "Spectral")),           
               prevalences = prevalences,
           detections = detections
           ) +
    labs(x = "Detection Threshold") 

Warning: The plot_core function is typically used with compositional 
                data. The data is not compositional. Make sure that you
                intend to operate on non-compositional data.

print(p)  

Subset entire relative abundance dataset with only core taxa

##identify core taxa
core.taxa <- core_members(pseq.rel, detection = 0, prevalence = 90/100)
ps.red.rel<-microbiome::transform(ps.red, "compositional")
ps.core<-subset_taxa(ps.red.rel, rownames(tax_table(ps.red.rel)) %in% core.taxa) #subsets relative abundance table from the ps.red sample dataset 
ps.core

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 51 taxa and 54 samples ]
sample_data() Sample Data:       [ 54 samples by 8 sample variables ]
tax_table()   Taxonomy Table:    [ 51 taxa by 7 taxonomic ranks ]
phy_tree()    Phylogenetic Tree: [ 51 tips and 50 internal nodes ]
refseq()      DNAStringSet:      [ 51 reference sequences ]

51 Taxa make up the Aiptasia “core” as defined by prevalence in 90 % of all control samples across time points at any relative abundance

range(sample_sums(ps.core))

[1] 0.4119657 0.9966101

range(sample_sums(subset_samples(ps.core, Condition=='Experimental')))

[1] 0.4119657 0.7929611

range(sample_sums(subset_samples(ps.core, Condition=='ABS1')))

[1] 0.7805256 0.9661604

range(sample_sums(subset_samples(ps.core, Condition=='ABS2')))

[1] 0.7163705 0.9966101

microbiome::plot_composition(ps.core, 
    plot.type = "barplot", sample.sort = "neatmap",group_by = 'Condition') 

core.taxa.df<-data.frame(tax_table(ps.core)) %>% arrange(Phylum, Class)
core.taxa.df

core.taxa.df %>% group_by(Class) %>% summarize(nr=n())

Fig 2B & C - betadiversity

# ordinate
set.seed(99)
ord<-ordinate(ps.core,"PCoA",distance = "bray")
ord.df<-plot_ordination(ps.core,ord,justDF=T)

fb<-ggplot(ord.df, aes(x=Axis.1, y =Axis.2, color=Condition, shape = Time, Fill=Condition)) + 
  geom_point(size=4) + 
  scale_shape_manual(values=shpsc) + 
  scale_color_manual(values=l2p16) + scale_fill_manual(values=l2p16) +
  theme_new #+ guides(shape=F, color=F)
fb

betadispersion

#run adonis
iDist <- phyloseq::distance(ps.core, method = "bray")  # note, must call phyloseq specifically if deseq2 is also loaded
pn_df <- data.frame(sample_data(ps.core))

#betadispersion
samps<-data.frame(sample_data(ps.red))

#betdisper
bd<-betadisper(iDist, paste(samps$Time,samps$Condition))
bd.df<-data.frame(distances=bd$distances, group1=bd$group,group2=bd$group)
bd.df<-bd.df %>% separate(group1,into = c('Time','Treatment'))
bd.df$Treatment<-factor(bd.df$Treatment, ordered=T, levels=c('Control','Experimental','ABS1','ABS2'))

c<-ggplot(bd.df, aes(group=group2, y=distances, x=Treatment, fill=Treatment)) + 
  geom_boxplot(alpha=0.5, outlier.shape = NA,color='black') + 
  scale_fill_manual(values=l2p16)+ scale_color_manual(values=l2p16)+
  geom_jitter(width=0.1,height=0,size=1, alpha=0.5) + facet_grid(~Time,scales='free',space='free')+ 
  theme_new + ylab('Distance to centriod') + xlab('') + guides(fill=F, group=F, color=F) + 
  theme(axis.text.x = element_blank(),panel.grid = element_blank(),plot.background = element_blank(),panel.grid.minor = element_blank(),rect  = element_blank())
c

stats

# Adonis test
mod<-adonis2(iDist ~ Condition*Time, data = pn_df)
mod

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = iDist ~ Condition * Time, data = pn_df)
               Df SumOfSqs      R2       F Pr(>F)    
Condition       3   2.5726 0.18752  7.7992  0.001 ***
Time            4   4.5488 0.33156 10.3425  0.001 ***
Condition:Time  4   1.9799 0.14431  4.5016  0.001 ***
Residual       42   4.6181 0.33661                   
Total          53  13.7194 1.00000                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

plot(mod)

#betadispersion
aov<-anova(bd)
aov

Analysis of Variance Table

Response: Distances
          Df  Sum Sq  Mean Sq F value Pr(>F)
Groups    11 0.44793 0.040721  1.2764 0.2711
Residuals 42 1.33996 0.031904               

bd

    Homogeneity of multivariate dispersions

Call: betadisper(d = iDist, group = paste(samps$Time, samps$Condition))

No. of Positive Eigenvalues: 36
No. of Negative Eigenvalues: 17

Average distance to median:
     T0 Control      T1 Control T1 Experimental         T4 ABS1         T4 ABS2      T4 Control         T5 ABS1         T5 ABS2 
         0.1688          0.1995          0.2540          0.1112          0.1887          0.2997          0.4210          0.3990 
     T5 Control         T6 ABS1         T6 ABS2      T6 Control 
         0.2133          0.1486          0.2152          0.2928 

Eigenvalues for PCoA axes:
(Showing 8 of 53 eigenvalues)
 PCoA1  PCoA2  PCoA3  PCoA4  PCoA5  PCoA6  PCoA7  PCoA8 
4.6413 2.5189 1.6697 0.9451 0.7118 0.5287 0.4723 0.3704 

tukey.bd<-TukeyHSD(bd, ordered=T)
tukey.bd

  Tukey multiple comparisons of means
    95% family-wise confidence level
    factor levels have been ordered

Fit: aov(formula = distances ~ group, data = df)

$group
                                  diff        lwr       upr     p adj
T6 ABS1-T4 ABS1            0.037406752 -0.4136519 0.4884654 1.0000000
T0 Control-T4 ABS1         0.057603553 -0.4466953 0.5619025 0.9999997
T4 ABS2-T4 ABS1            0.077534137 -0.3735245 0.5285928 0.9999773
T1 Control-T4 ABS1         0.088337882 -0.3833905 0.5600663 0.9999466
T5 Control-T4 ABS1         0.102137167 -0.3695913 0.5738656 0.9997794
T6 ABS2-T4 ABS1            0.104036920 -0.3470217 0.5550956 0.9995975
T1 Experimental-T4 ABS1    0.142810858 -0.3082478 0.5938695 0.9933197
T6 Control-T4 ABS1         0.181585017 -0.2694736 0.6326437 0.9587465
T4 Control-T4 ABS1         0.188454688 -0.2626040 0.6395133 0.9469959
T5 ABS2-T4 ABS1            0.287766418 -0.1632922 0.7388251 0.5545898
T5 ABS1-T4 ABS1            0.309798021 -0.1412606 0.7608567 0.4439119
T0 Control-T6 ABS1         0.020196802 -0.4308618 0.4712554 1.0000000
T4 ABS2-T6 ABS1            0.040127385 -0.3505009 0.4307556 0.9999999
T1 Control-T6 ABS1         0.050931130 -0.3633927 0.4652550 0.9999993
T5 Control-T6 ABS1         0.064730416 -0.3495934 0.4790542 0.9999914
T6 ABS2-T6 ABS1            0.066630169 -0.3239981 0.4572584 0.9999790
T1 Experimental-T6 ABS1    0.105404106 -0.2852241 0.4960324 0.9983034
T6 Control-T6 ABS1         0.144178266 -0.2464500 0.5348065 0.9778134
T4 Control-T6 ABS1         0.151047936 -0.2395803 0.5416762 0.9688830
T5 ABS2-T6 ABS1            0.250359666 -0.1402686 0.6409879 0.5478056
T5 ABS1-T6 ABS1            0.272391269 -0.1182370 0.6630195 0.4211407
T4 ABS2-T0 Control         0.019930583 -0.4311281 0.4709892 1.0000000
T1 Control-T0 Control      0.030734329 -0.4409941 0.5024628 1.0000000
T5 Control-T0 Control      0.044533614 -0.4271948 0.5162620 1.0000000
T6 ABS2-T0 Control         0.046433367 -0.4046253 0.4974920 0.9999999
T1 Experimental-T0 Control 0.085207305 -0.3658513 0.5362660 0.9999417
T6 Control-T0 Control      0.123981464 -0.3270772 0.5750401 0.9980013
T4 Control-T0 Control      0.130851135 -0.3202075 0.5819098 0.9968010
T5 ABS2-T0 Control         0.230162864 -0.2208958 0.6812215 0.8272338
T5 ABS1-T0 Control         0.252194468 -0.1988642 0.7032531 0.7324375
T1 Control-T4 ABS2         0.010803745 -0.4035201 0.4251276 1.0000000
T5 Control-T4 ABS2         0.024603031 -0.3897208 0.4389269 1.0000000
T6 ABS2-T4 ABS2            0.026502784 -0.3641255 0.4171310 1.0000000
T1 Experimental-T4 ABS2    0.065276721 -0.3253515 0.4559050 0.9999830
T6 Control-T4 ABS2         0.104050881 -0.2865774 0.4946791 0.9984881
T4 Control-T4 ABS2         0.110920551 -0.2797077 0.5015488 0.9973420
T5 ABS2-T4 ABS2            0.210232281 -0.1803960 0.6008605 0.7752301
T5 ABS1-T4 ABS2            0.232263884 -0.1583644 0.6228921 0.6544875
T5 Control-T1 Control      0.013799285 -0.4229364 0.4505349 1.0000000
T6 ABS2-T1 Control         0.015699038 -0.3986248 0.4300229 1.0000000
T1 Experimental-T1 Control 0.054472976 -0.3598508 0.4687968 0.9999986
T6 Control-T1 Control      0.093247135 -0.3210767 0.5075710 0.9996809
T4 Control-T1 Control      0.100116806 -0.3142070 0.5144406 0.9993781
T5 ABS2-T1 Control         0.199428536 -0.2148953 0.6137524 0.8738175
T5 ABS1-T1 Control         0.221460139 -0.1928637 0.6357840 0.7824422
T6 ABS2-T5 Control         0.001899753 -0.4124241 0.4162236 1.0000000
T1 Experimental-T5 Control 0.040673691 -0.3736501 0.4549975 0.9999999
T6 Control-T5 Control      0.079447850 -0.3348760 0.4937717 0.9999324
T4 Control-T5 Control      0.086317521 -0.3280063 0.5006413 0.9998478
T5 ABS2-T5 Control         0.185629251 -0.2286946 0.5999531 0.9172050
T5 ABS1-T5 Control         0.207660854 -0.2066630 0.6219847 0.8426829
T1 Experimental-T6 ABS2    0.038773938 -0.3518543 0.4294022 0.9999999
T6 Control-T6 ABS2         0.077548097 -0.3130802 0.4681763 0.9999049
T4 Control-T6 ABS2         0.084417768 -0.3062105 0.4750460 0.9997834
T5 ABS2-T6 ABS2            0.183729497 -0.2068988 0.5743577 0.8893728
T5 ABS1-T6 ABS2            0.205761101 -0.1848671 0.5963893 0.7973195
T6 Control-T1 Experimental 0.038774159 -0.3518541 0.4294024 0.9999999
T4 Control-T1 Experimental 0.045643830 -0.3449844 0.4362721 0.9999996
T5 ABS2-T1 Experimental    0.144955560 -0.2456727 0.5355838 0.9769141
T5 ABS1-T1 Experimental    0.166987163 -0.2236411 0.5576154 0.9384202
T4 Control-T6 Control      0.006869671 -0.3837586 0.3974979 1.0000000
T5 ABS2-T6 Control         0.106181400 -0.2844468 0.4968096 0.9981889
T5 ABS1-T6 Control         0.128213004 -0.2624152 0.5188413 0.9910467
T5 ABS2-T4 Control         0.099311730 -0.2913165 0.4899400 0.9990076
T5 ABS1-T4 Control         0.121343333 -0.2692849 0.5119716 0.9942930
T5 ABS1-T5 ABS2            0.022031603 -0.3685966 0.4126599 1.0000000

#replace NAs with 0s to get multcomp to run 
tukey$`Condition:Time`[!complete.cases(tukey$`Condition:Time`),] <- 0

library(multcompView)
tukey.cld <-multcompLetters4(aov,comp = tukey.bd)
print(tukey.cld)
multcompLetters(tukey.bd$group[,4])

#Fig 2D plot distance tree

tr<-plot_tree(ps.core, 'treeonly',label.tips="taxa_names",nodelabf=nodeplotblank,ladderize="left", justify = "yes please")

Warning: Some columns are a multi-column type (such as a matrix column): [1, 2, 3, 4, 5, 6, 7]. setDT will retain these columns as-is but subsequent operations like grouping and joining may fail. Please consider as.data.table() instead which will create a new column for each embedded column.Warning: Some columns are a multi-column type (such as a matrix column): [1, 2, 3, 4, 5, 6, 7]. setDT will retain these columns as-is but subsequent operations like grouping and joining may fail. Please consider as.data.table() instead which will create a new column for each embedded column.

asv.order<-data.frame(tr$data) %>% arrange(y) 
asv.order<-na.omit(as.vector(asv.order$OTU))
asv.order

 [1] "ASV1"   "ASV5"   "ASV56"  "ASV3"   "ASV19"  "ASV6"   "ASV86"  "ASV4"   "ASV9"   "ASV25"  "ASV24"  "ASV21"  "ASV61"  "ASV110"
[15] "ASV44"  "ASV22"  "ASV35"  "ASV39"  "ASV16"  "ASV37"  "ASV67"  "ASV13"  "ASV152" "ASV51"  "ASV95"  "ASV15"  "ASV81"  "ASV27" 
[29] "ASV47"  "ASV34"  "ASV8"   "ASV90"  "ASV75"  "ASV26"  "ASV12"  "ASV14"  "ASV32"  "ASV79"  "ASV87"  "ASV2"   "ASV49"  "ASV91" 
[43] "ASV40"  "ASV101" "ASV43"  "ASV38"  "ASV59"  "ASV7"   "ASV20"  "ASV17"  "ASV33" 
attr(,"na.action")
 [1]  4  6  8 10 11 13 17 19 20 21 23 25 27 31 32 34 36 37 39 40 44 46 48 50 52 54 55 56 58 60 63 65 67 70 71 73 76 77 79 81 84 85
[43] 87 89 91 93 95 97 99
attr(,"class")
[1] "omit"

d<-ggtree::ggtree(ps.core) + ggtree::geom_tiplab(align=T, hjust=-.3) + ggtree::hexpand(.01) 

! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.

d

Error in `geom_segment2()`:
! Problem while converting geom to grob.
ℹ Error occurred in the 3rd layer.
Caused by error in `check.length()`:
! 'gpar' element 'lwd' must not be length 0
Backtrace:
  1. base (local) `<fn>`(x)
  2. ggplot2:::print.ggplot(x)
  4. ggplot2:::ggplot_gtable.ggplot_built(data)
  5. ggplot2:::by_layer(...)
 12. ggplot2 (local) f(l = layers[[i]], d = data[[i]])
     ...
 24. grid::gpar(...)
 25. grid:::validGP(list(...))
 26. grid (local) numnotnull("lwd")
 27. grid (local) check.length(gparname)
 28. base::stop(...)

library(ggtree)

ggtree v3.4.4 For help: https://yulab-smu.top/treedata-book/

If you use the ggtree package suite in published research, please cite the appropriate paper(s):

Guangchuang Yu, David Smith, Huachen Zhu, Yi Guan, Tommy Tsan-Yuk Lam. ggtree: an R package for visualization and
annotation of phylogenetic trees with their covariates and other associated data. Methods in Ecology and Evolution.
2017, 8(1):28-36. doi:10.1111/2041-210X.12628

LG Wang, TTY Lam, S Xu, Z Dai, L Zhou, T Feng, P Guo, CW Dunn, BR Jones, T Bradley, H Zhu, Y Guan, Y Jiang, G Yu.
treeio: an R package for phylogenetic tree input and output with richly annotated and associated data. Molecular
Biology and Evolution. 2020, 37(2):599-603. doi: 10.1093/molbev/msz240

Guangchuang Yu, Tommy Tsan-Yuk Lam, Huachen Zhu, Yi Guan. Two methods for mapping and visualizing associated data on
phylogeny using ggtree. Molecular Biology and Evolution. 2018, 35(12):3041-3043. doi:10.1093/molbev/msy194 

Attaching package: ‘ggtree’

The following object is masked from ‘package:Biostrings’:

    collapse

The following object is masked from ‘package:IRanges’:

    collapse

The following object is masked from ‘package:S4Vectors’:

    expand

The following object is masked from ‘package:ape’:

    rotate

The following object is masked from ‘package:tidyr’:

    expand

ggtree(ps.core) + geom_tiplab(align=F, hjust=-.3, aes(label=Class))

! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.

Fig 2E

# Bubble plot of relative abundances in ps.core
# Include only 1 control sample per time point to reduce complexity of plot 
ps.core.control<-subset_samples(ps.core,Condition=='Control') %>%  subset_samples(Sample_names %in%c('SAM_133_S1','SAM_157_S10','SAM_347_S56','SAM_364_S73','SAM_409_S89'))

Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’

ps.core.other<-subset_samples(ps.core,!Condition=='Control') 

Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’

#merge
ps.core.new<-merge_phyloseq(ps.core.control,ps.core.other)

Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’

# place samples into correct order
control.samps<-data.frame(sample_data(ps.core.control))$Sample_names
other.samps<-data.frame(sample_data(ps.core.other))$Sample_names
samps.order<-c(control.samps,other.samps)

#Add Grouping variable for plot 
samps.new<-data.frame(sample_data(ps.core.new)) %>% mutate(PlotGroup=case_when(Condition=='Control' ~ 'ASW', 
                                                                               Condition=='Experimental' ~ 'FASW 33 d',
                                                                               Condition=='ABS1' & Time=='T4' ~ "ABS1 55 d",
                                                                                Condition=='ABS1' & Time=='T5' ~ "ABS1 61 d",
                                                                                Condition=='ABS1' & Time=='T6' ~ "ABS1 76 d",
                                                                                Condition=='ABS2' & Time=='T4' ~ "ABS2 55 d",
                                                                                Condition=='ABS2' & Time=='T5' ~ "ABS2 61 d",
                                                                                Condition=='ABS2' & Time=='T6' ~ "ABS2 76 d"))

sample_data(ps.core.new)<-samps.new

Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’

#make color palette for core taxa via family level classification 
df.taxa<-data.frame(ASV=rownames(tax_table(ps.core.new)),tax_table(ps.core.new)) %>% arrange(Phylum, Class, Order, Family, Genus)
df.taxa<-data.frame(ASV=rownames(tax_table(ps.core.new)),tax_table(ps.core.new)) %>% arrange(factor(ASV, levels=asv.order))
df.taxa$Family2<-ifelse(is.na(df.taxa$Family),paste(df.taxa$Class,'-unknown',sep=''),df.taxa$Family) #fix na family level classificaiton 
df.taxa$Family2<-ifelse(df.taxa$Family2=='NA-unknown',paste(df.taxa$Phylum,'-unknown',sep=''),df.taxa$Family2) #fix na family level classificaiton 
df.taxa

#set color pallet based on taxa order in tree, reds gammas, yellows, alphas, purples 
clr<-c('#FCD9DF','#38050D','#F7A1AF','#F05670','#EB1E40','#CE1231','#960D24', #"Alteromonadaceae","Gammaproteobacteria-unknown","Vibrionaceae","Halomonadaceae","Pseudomonadaceae","Oleiphilaceae",Spongiibacteraceae 
  '#D46021',# Akkermansiaceae
  'mediumpurple4','mediumpurple1', #"Phycisphaerae-unknown","Pirellulaceae" 
   '#9BC0AF', #Nannocystaceae
  '#187B4D', #Mycobacteriaceae
  '#093149','#49AEE9','#A4D7F4',#"Bacteroidia-unknown","Saprospiraceae","Flavobacteriaceae" 
  '#F9F5DC',"#F0E7A8",'#F7EB8D',"#E7DA73",#"Emcibacteraceae","Hyphomonadaceae","GCA-2696645","Sphingomonadaceae" 
  "gray",#"Proteobacteria-unknown"
  '#F3E04B','#574F0F','#F0DA32','#E4CB11','#BEA90E','#85760A'#"Devosiaceae","Stappiaceae","Alphaproteobacteria-unknown", "Cohaesibacteraceae"
  #"Rhodobacteraceae","Rhizobiaceae" 
)    


names(clr)<-unique(df.taxa$Family2)

make bubble plot

# get relavent dataframes 
md<-data.frame(sample_data(ps.core.new))
df.asv<-data.frame(Sample_names=rownames(otu_table(ps.core.new)),otu_table(ps.core.new))

#make long dataframe for plotting
df.asv.long<-pivot_longer(df.asv,cols=!Sample_names,names_to='ASV',values_to = 'RelAbund') #melt dataframe
df.relabund<-df.asv.long %>% left_join(y=md, by='Sample_names') #add in metadata
df.relabund<-df.relabund %>% left_join(y=df.taxa, by='ASV') #add in taxa info
df.relabund<-df.relabund[df.relabund$RelAbund>0,]

#re-order ASVs
df.relabund$ASV<-factor(df.relabund$ASV, ordered=T, levels=asv.order)

#reoder samples 
df.relabund$Sample_names<-factor(df.relabund$Sample_names, ordered=T, levels=samps.order)

#re-order treatments
unique(df.relabund$PlotGroup)

[1] "ASW"       "FASW 33 d" "ABS1 55 d" "ABS2 55 d" "ABS1 61 d" "ABS2 61 d" "ABS1 76 d" "ABS2 76 d"

df.relabund$PlotGroup<-factor(df.relabund$PlotGroup, ordered=T, levels=c("ASW","FASW 33 d", "ABS1 55 d", "ABS2 55 d", "ABS1 61 d", "ABS2 61 d", "ABS1 76 d","ABS2 76 d"))

head(df.relabund)

e<-ggplot(df.relabund, aes(x = Sample_names, y = ASV, fill = Family2)) +
  geom_point(aes(size = RelAbund*10), alpha = 0.8, shape = 21) + 
  scale_size_continuous(limits = c(0, 70), range = c(0,20), breaks = c(0.01, 5, 10, 20, 40, 80)) +
  theme_new + xlab('') +ylab('') + 
  facet_grid(~PlotGroup, scales='free',space='free')+
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
  labs(size = "Relative Abundance", color = "Family")  + 
  scale_fill_manual(values=clr) + guides(fill=F)
e

#ggsave('maggie_plots/bubbles-with-legend.pdf')

put all together

library(patchwork)

layout3<-"AABBBB
CCBBBB
DEEEEE
DEEEEE
DEEEEE
"

e2<-e+theme(panel.grid = element_blank(),rect = element_blank(),panel.grid.minor = element_blank(),panel.grid.major = element_blank())
e2

fb2<-fb + guides(color="none", fill="none", shape="none")

a + fb2 + c + d + e2 + plot_layout(design=layout3)

ggsave('Figures/Figure2-new.pdf', width=11,height=10)

#save core taxa as .csv
df.taxa<-df.taxa %>% unite('string',c(Phylum,Class,Order,Family,Genus,Species), sep='_;',na.rm=F, remove=F)

head(df.taxa)

write.csv(df.taxa, 'Figures/core.taxa.csv')

Figure S1

samps<-data.frame(sample_data(ps.core))
samps$Condition<- factor(samps$Condition, ordered=T, levels=c('Control', 'Experimental','ABS1','ABS2'))
samps<-samps %>% arrange(Condition,Time)

#make new grouping factor 
samps$grp<-paste(samps$Condition, samps$Time)
samps$grp<-factor(samps$grp, ordered = T, levels=c("Control T0","Control T1","Control T4", "Control T5","Control T6","Experimental T1", "ABS1 T4","ABS1 T5","ABS1 T6","ABS2 T4", "ABS2 T5","ABS2 T6"))
sample_data(ps.core)<-samps

ps.core %>%
  aggregate_taxa(level='Family') %>% 
  plot_composition(plot.type='barplot',group_by = 'grp',sample.sort = 'neatmap')  +
  scale_fill_manual(values=clr)+
  scale_y_continuous(label = scales::percent) + 
  theme_new + guides(fill=F)+
  # Removes sample names and ticks
  theme( axis.text.x = element_blank(), axis.ticks.x=element_blank(), plot.background = element_blank()) 
ggsave('Figures/figS1.pdf',height=6,width=12)